2 November 2016

Urban Observatory Air Quality Data

Urban Observatory Readings

Bayesian Parametric Inference

  • Given data, \(y\), we want to fit a model with parameters \(\theta\)

  • Bayes theorem allows us to do this

\[p(\theta | y) = \frac{p(\theta)p(y | \theta)}{p(y)} \]

  • \(p(\theta | y)\) is the posterior distribution, \(p(\theta)\) is the prior distribution, \(p(y|\theta\)) the likelihood and \(p(y) = \int_\theta p(y | \theta) p(\theta) \textrm{d}\theta\)

  • Once parameters have been determined, the parametric model can be used to answer questions about the process which generated the data, \(y\)

Streaming Data

  • Streaming data is observed as a chronological sequence of values, \(Y\), with an associated timestamp, \(t_i\)

\[ Y(t_{0:N}) = \{y(t_0), y(t_1), \dots , y(t_N) \} \]

  • Typically streaming data is unbounded:

\[y(t_0), y(t_1), \dots \]

  • We want to forecast future observations, which have not yet been observed

  • We do this by fitting a parametric model which describes the process

Irregularly Observed Data

  • Streaming data is often irregularly observed, consider first a completely observed discrete sensor for a lightbulb

  • The lightbulb can be considered either on, or off, here is a regularly observed series of the lightbulb

  • The lightbulb reports its status every five minutes, leading to \(\frac{60}{5} \times 24 = 288\) observations each day

Discrete Irregularly Observed Data

  • Recording the lightbulb data irregularly can increase accuracy, and reduce storage costs

  • In the regularly observed data there are 288 observations a day, in the irregular case there is one observation

  • With irregularly observed data, we can accurately answer questions like: "at 12pm on Tuesday, how many lights were on in The Core?"

Modelling Irregularly Observed Data

POMP Models

  • Partially Observed Markov Process model

\[\begin{align*} y(t_i)|\eta(t_i) &\sim \pi(y(t_i) | \eta(t_i)), \\ \eta(t_i)|\textbf{x}(t_i) &= g(F_{t_i}^T \textbf{x}(t_i)), \\ \textbf{X}(t_i) | \textbf{x}(t_{i-1}) &\sim p(\textbf{x}(t_i) | \textbf{x}(t_{i-1})) \end{align*}\]

  • The observation distribution \(\pi(.)\) is flexible
  • The function \(g\) permits a deterministic non-linear transformation of the state space
  • \(F_t\) is a time dependent vector representing a linear transformation
  • The state space, \(\textbf{x}(t)\) is a continuous time Markov process

POMP Model

Representation of a POMP model as a directed acyclic graph (DAG)

The State Space

Diffusion Process Plot, Generalised Brownian Motion

The State Space

Diffusion Process Plot, Ornstein-Uhlenbeck Process

Example POMP Model

A Poisson Model

\[\begin{align*} N(t) &\sim \textrm{Poisson}(N(t) | \lambda(t)) \\ \lambda(t_i) &= \exp\{ x(t) \} \\ \textrm{d}X(t) &= \mu \textrm{d}t + \sigma \textrm{d}W(t) \end{align*}\]

Composing POMP models

Modelling Complex Processes

  • Define an associative binary composition operator for models
  • The Left-hand observation model and linking function, \(g\) are chosen
  • The Linear transformation vectors are concatenated: \(F^{(3)}_t = \begin{bmatrix} F^{(1)}_t \\ F^{(2)}_t \end{bmatrix}\)
  • The latent state of each model advances according to its state evolution equation:

\[ \textbf{X}(t_i) | \textbf{x}(t_{i-1}) \sim \begin{pmatrix} p_1(\textbf{x}^{(1)}(t_i) | \textbf{x}^{(1)}(t_{i-1})) \\ p_2(\textbf{x}^{(2)}(t_i) | \textbf{x}^{(2)}(t_{i-1})) \end{pmatrix} \]

  • The class of unparameterised models and the composition operator forms a semi group

Composed POMP Example

A Seasonal Poisson Model

  • First define a seasonal model, with an Ornstein-Uhlenbeck state space

\[\begin{align*} y(t) &\sim \mathcal{N}(y(t) | \mu(t), \sigma) \\ \mu(t) &= F_t^T \textbf{x}(t) \\ \textrm{d}\textbf{X}(t) &= \alpha(\theta - \textbf{x}(t)) \textrm{d}t + \Sigma \textrm{d}W(t) \end{align*}\]

  • \(F_t^T\) is a time dependent vector of fourier components representing seasonality:

\[F_t = \begin{pmatrix} \cos(\omega t) \\ \sin(\omega t) \\ \cos(2\omega t) \\ \sin(2\omega t) \\ \cos(3\omega t) \\ \sin(3\omega t) \end{pmatrix}\]

Composed POMP Example

  • Define a Poisson model as previously:

\[\begin{align*} N(t) &\sim \textrm{Poisson}(N(t) | \lambda(t)) \\ \lambda(t_i) &= \exp\{ x(t) \} \\ \textrm{d}X(t) &= \mu \textrm{d}t + \sigma \textrm{d}W(t) \end{align*}\]

Composed Pomp Example

A Seasonal-Poisson Model

\[\begin{align*} N(t) &\sim \textrm{Poisson}(N(t) | \lambda(t)) \\ \lambda(t_i) &= \exp\{ F_t^T x(t) \} \\ \textbf{X}(t_i)|\textbf{x}(t_{i-1}) &\sim \begin{pmatrix} p_1(\textbf{x}^{(1)}(t_i) | \textbf{x}^{(1)}(t_{i-1})) \\ p_2(\textbf{x}^{(2)}(t_i) | \textbf{x}^{(2)}(t_{i-1})) \end{pmatrix} \end{align*}\]

Programming with Streams

Functional Streams

  • Streams can provide a way of working with unbounded data using finite resources
  • Functional streams are modular, composable and easy to reason about

Stream

An infinite list

  • Streams can be thought of as an infinite list
  • A Stream is pair, with the left hand value being the first element of the stream and the right hand value a thunk (a function with no arguments)
  • Define a (infinite) stream of natural numbers using from
scala> val naturalNumbers = Stream.from(1)

Stream(1, ?)

Akka Streams

  • Akka Streams: A Scala library for stream processing

  • Akka streams have three main abstractions:

  • A Source is a definition of a stream, it can be a "pure" stream, or a database or webservice call

  • A Flow is a processing stage, which can be used to transform a Source to another Source

  • A Sink defines what happens at the end of the stream, usually an effect such as writing to a file

Higher Order Functions

Operations on Streams

  • foldLeft (or foldRight) can be used (with a seed) to reduce a foldable data structure by recursively applying a binary operation
scala> Stream.from(1).
take(10).
foldLeft(0)(_ + _)

55

  • foldLeft will reduce from the left, ((1 + 2) + 3) + 4

  • foldRight will reduce from the right, 1 + (2 + (3 + 4))

  • reduce is equivalent to foldLeft for associative and commutative binary operations and can be applied in parallel

Higher Order Functions

Operations on Streams

  • scanLeft can be used to accumulate a running sum:
scala> Stream.from(1).
take(10).
scanLeft(0)(_ + _).
foreach(n => print(s"$n, "))

0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55

POMP models as a Stream

  • The latent state is a Markov process, which means future observations only depend on the current observation

\[p(x(t_n) | x(t_{n-1}), \dots x(t_0)) = p(x(t_n) | x(t_{n-1}))\]

  • In order to simulate the Markov Process, we can use the dual of fold, unfold, to simulate a random walk:
val x0 = Gaussian(0.0, 1.0).draw

Source.unfold(x0)(a => Some((Gaussian(a, sigma).draw, a)))

Bayesian Inference with Streams

The Bootstrap Particle Filter

  • In general, POMP models require simulation based filtering (such as the bootstrap filter below) to determine the behaviour of the latent state:
  1. Assume we have a collection of weighted "particles" representing the latent state at time \(t_{k-1}\) \(\{x^{(i)}(t_{k-1}), w^{(i)}(t_{k-1}), i = 1,\dots, N\}\)
  2. Advance the particles to the time of the next observation, using the state transition model \(x^{(i)}(t) = p(x^{(i)}(t_k) | x^{(i)}(t_{k-1}))\)
  3. Calculate the likelihood (weight) of each particle, given the observation at time \(t_k\), \(w^{(i)}(t_k) = p(y(t_k) | x(t_k))\)
  4. Resample the particles according to the weights, ie. select particle \(j\) with probability \(w^{(j)}(t_k)\)
  5. Add the resampled state to the state vector \(x(t_{1:k}) = (x(t_{1:k-1}), x(t_{k}))\), set \(i = i + 1\) and go to 2

Illustration of Bootstrap Particle Filter

  • The solid line represents the simulated state, the points represent the particle cloud with the opacity representing each particles likelihood

Bayesian Inference With Streams

  • In order to perform the bootstrap particle filter, we need to know the time of the previous observation and the previous particle cloud

  • First define the observation datatype and the state of the bootstrap particle filter

    case class Data(time: Datetime, observation: Double)
    case class FilterState(t0: Datetime, state: List[State])
  • The signature of the filtering step function is: scala val filterStep: (Data, FilterState) => FilterState
  • The the filter can by ran using an akka Flow and the scan operation, the initial state init contains the particle cloud and the time at the start of the application of the filter

    def filter(init: FilterState) = Flow[Data].scan(init)(filterStep)

Conclusion

  • Bayesian Inference for irregular time series
  • Functional Streams help, by handling unbounded data in memory

  • I have wrote a scala library which implements composable POMP models and inference methods
  • Look at the code: git.io/statespace
  • Read the paper: arXiv:1609.00635